This paper aims to explore changes in the tempo of cultural evolutionary changes using range data (richness with start and end dates). The decision to focus on range data is to highlight that this type of data is very common in many forms of evolutionary analysis, as taphonomic processes inherently reduce the granularity of data (this is the focus of the collaborative paper). This not only includes paleobiology, but also archaeology and anthropology, and a range of other sciences as well.
Here, we suggest that the most effective set of tools for handling range data, is a model-based approach (like that implemented in PyRate), as opposed to relying on empirical rates or raw richness counts. The reasons for this are two-fold: 1) Similar patterns in richness can be the result of different underlying processes and 2) empirical diversification rates are often very noisy, and it is difficult to tell whether changes in empirical rates over time are signficant. Based on this, the key questions of the paper would be:
The strucutre of the analysis include the following:
Calculate the expected changes in richness and diversification rates (origination, extinction, net diversification, turnover) based on the expectations of the netural model
Calculate the empirical values and confidence intervals of richness and diversification rates from a large number of replicated datasets
Estimate diversification rates using an underlying birth-death model of diversification (i.e PyRate).
Compare the expected, empirical, and modeled results. The expected values require us to know the frequencies of different cultural variants, whereas the empirical and modeled results only require range data (richness plus dates)
TBD
TBD
library(tidyverse)
source("~/Dropbox/cult_tempo/cult_tempo_functions.R")
#source from GitHub once committed
#Global model settings
ts = 50 # number of timesteps
warm = 500 #number of warmup timesteps
rep = 500 #number of neutral model replicates to create
#Population Size
N1 = 1000 #starting population size
N2 = 3000 #ending population size
#Innovation Rate
mu1 = 0.01 #starting innovation rate
mu2 = 0.01 #ending innovation rate
#Strength of Frequency Bias Transmission (Under Development)
b1 = 0 #starting strength of transmission
b2 = 0 #ending strength of transmission
pop_sizes <- tibble(timesteps = seq(ts, 1, -1)) %>%
mutate(linear=seq(from = N1, to = N2, length.out = length(timesteps))) %>%
mutate(logistic = sig_growth(timesteps, assympote = N2, lower_bound = N1, midpoint = median(timesteps), scale = 0.2)) %>%
mutate(exponetial = expo_growth(timesteps, upper_bound = N2, lower_bound = N1, scale = 0.1)) %>%
mutate(cycle = cyclical(timesteps, midpoint = (N2+N1)/2, amplitude = N2-(N2+N1)/2, scale = 4))
#Population Size
N1 = 1000 #starting population size
N2 = 1000 #ending population size
#Innovation Rate
mu1 = 0.01 #starting innovation rate
mu2 = 0.05 #ending innovation rate
mu_shifts <- seq(from=1,to=ts,by=1)
mu_values <- seq(from=mu1, to=mu2, length.out = ts)
innov_linear <- innovation_rate(mu_shifts = mu_shifts, mu_values=mu_values, ts = ts)
mu_shifts <- c(1,25)
mu_values <- c(0.01,0.05)
innov_key_invention <- innovation_rate(mu_shifts = mu_shifts, mu_values=mu_values, ts = ts)
mu_shifts <- c(1,11,21,31,41)
mu_values <- c(0.01, 0.02, 0.03, 0.04, 0.05)
innov_incremental <- innovation_rate(mu_shifts = mu_shifts, mu_values=mu_values, ts = ts)
innov_rates <- tibble(timesteps = seq(ts, 1, -1)) %>%
mutate(linear = innov_linear) %>%
mutate(invention = innov_key_invention) %>%
mutate(incremental = innov_incremental)
#Strength of Frequency Bias Transmission (Under Development)
b1 = 0 #starting strength of transmission
b2 = 0 #ending strength of transmission
#Population Size
N1 = 1000 #starting population size
N2 = 1000 #ending population size
#Innovation Rate
mu1 = 0.01 #starting innovation rate
mu2 = 0.01 #ending innovation rate
#Strength of Frequency Bias Transmission (Under Development)
b1 = -0.03 #starting strength of transmission
b2 = 0.03 #ending strength of transmission
cult_shifts <- seq(from=1,to=ts,by=1)
b_values <- seq(from=b1, to=b2, length.out = ts)
trans_linear <- trans_rate(cult_shifts = cult_shifts, b_values = b_values, ts = ts)
cult_shifts <- c(1,25)
b_values <- c(-0.03,0.03)
trans_single <- trans_rate(cult_shifts = cult_shifts, b_values=b_values, ts = ts)
cult_shifts <- c(1,11,21,31,41)
b_values <- c(-0.03, -0.01, 0, 0.01, 0.03)
trans_multiple <- trans_rate(cult_shifts = cult_shifts, b_values=b_values, ts = ts)
#One Replicate
neutral_linear <- neutral(N = pop_sizes$linear, mu = mu, warm = warm, ts = ts)
neutral_logistic <- neutral(N = pop_sizes$logistic, mu = mu, warm = warm, ts = ts)
neutral_exponential <- neutral(N = pop_sizes$exponetial, mu = mu, warm = warm, ts = ts)
neutral_cycle <- neutral(N = pop_sizes$cycle, mu = mu, warm = warm, ts = ts)
These functions make multiple replicates in order to identify the envelope of empirical values. These are also saved as a RData object and can be easily loaded in for plotting without having to go through the replicates procedure.
#Multiple Replicates
neutral_linear_reps <- neutral_replicates(N = pop_sizes$linear, mu = mu, warm = warm, ts = ts, reps = rep)
neutral_linear_empirical<-lapply(neutral_linear_reps,empirical_values)
neutral_linear_summary <- empirical_summary(neutral_linear_empirical)
neutral_logistic_reps <- neutral_replicates(N = pop_sizes$logistic, mu = mu, warm = warm, ts = ts, reps = rep)
neutral_logistic_empirical<-lapply(neutral_logistic_reps,empirical_values)
neutral_logistic_summary <- empirical_summary(neutral_logistic_empirical)
neutral_exponential_reps <- neutral_replicates(N = pop_sizes$exponetial, mu = mu, warm = warm, ts = ts, reps = rep)
neutral_exponential_empirical<-lapply(neutral_exponential_reps,empirical_values)
neutral_exponential_summary <- empirical_summary(neutral_exponential_empirical)
neutral_cycle_reps <- neutral_replicates(N = pop_sizes$cycle, mu = mu, warm = warm, ts = ts, reps = rep)
neutral_cycle_empirical<-lapply(neutral_cycle_reps,empirical_values)
neutral_cycle_summary <- empirical_summary(neutral_cycle_empirical)
save(ts,warm,rep,N,mu,b,
pop_sizes,
neutral_linear_summary,
neutral_logistic_summary,
neutral_exponential_summary,
neutral_cycle_summary,exp.std.div,
file="neutural_pop_data.RData")
N = rep(N1,ts)
neu_innov_linear <- neutral(N = N, mu = innov_rates$linear, warm = warm, ts = ts)
neu_innov_invention <- neutral(N = N, mu = innov_rates$invention, warm = warm, ts = ts)
neu_innov_incremental <- neutral(N = N, mu = innov_rates$incremental, warm = warm, ts = ts)
neu_innov_linear_reps <- neutral_replicates(N = N, mu = innov_rates$linear, warm = warm, ts = ts, reps = rep)
neu_innov_linear_empirical<-lapply(neu_innov_linear_reps,empirical_values)
neu_innov_linear_summary <- empirical_summary(neu_innov_linear_empirical)
neu_innov_invention_reps <- neutral_replicates(N = N, mu = innov_rates$invention, warm = warm, ts = ts, reps = rep)
neu_innov_invention_empirical<-lapply(neu_innov_invention_reps,empirical_values)
neu_innov_invention_summary <- empirical_summary(neu_innov_invention_empirical)
neu_innov_incremental_reps <- neutral_replicates(N = N, mu = innov_rates$incremental, warm = warm, ts = ts, reps = rep)
neu_innov_incremental_empirical<-lapply(neu_innov_incremental_reps,empirical_values)
neu_innov_incremental_summary <- empirical_summary(neu_innov_incremental_empirical)
save(ts,warm,rep,N,
innov_rates,
neu_innov_linear_summary,
neu_innov_invention_summary,
neu_innov_incremental_summary,
file="neutral_innov_data.RData")
TBD
#Population Growth
neutral_linear_LiteRate<-matrix_to_literate(neutral_linear)
write.table(neutral_linear_LiteRate,"neutral_linear_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
neutral_logistic_LiteRate<-matrix_to_literate(neutral_logistic)
write.table(neutral_logistic_LiteRate,"neutral_logistic_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
neutral_exponential_LiteRate<-matrix_to_literate(neutral_exponential)
write.table(neutral_exponential_LiteRate,"neutral_exponential_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
neutral_cycle_LiteRate<-matrix_to_literate(neutral_cycle)
write.table(neutral_cycle_LiteRate,"neutral_cycle_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
#Innovation Rates
neu_innov_linear_LiteRate<-matrix_to_literate(neu_innov_linear)
write.table(neu_innov_linear_LiteRate,"neu_innov_linear_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
neu_innov_invention_LiteRate<-matrix_to_literate(neu_innov_invention)
write.table(neu_innov_invention_LiteRate,"neu_innov_invention_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
neu_innov_incremental_LiteRate<-matrix_to_literate(neu_innov_incremental)
write.table(neu_innov_incremental_LiteRate,"neu_innov_incremental_LiteRate.txt",quote = FALSE, sep="\t",row.names = FALSE)
#Frequency Transmission - TBD
#Population Growth Model
#-n 10000000 (default number of mcmc iterations)
#-s 1000 (default sampling frequency)
#-TBP means dates are in Time Before Present)
python3 LiteRateForward.py -d ~/neutral_linear_LiteRate.txt -TBP #population linear
python3 LiteRateForward.py -d ~/neutral_logistic_LiteRate.txt -TBP #population linear
python3 LiteRateForward.py -d ~/neutral_exponential_LiteRate.txt -TBP #population linear
python3 LiteRateForward.py -d ~/neutral_cycle_LiteRate.txt -TBP #population linear
#Innovation Rate Models
python3 LiteRateForward.py -d ~/neu_innov_linear_LiteRate.txt -TBP #innovation linear
python3 LiteRateForward.py -d ~/neu_innov_invention_LiteRate.txt -TBP #innovation linear
python3 LiteRateForward.py -d ~/neu_innov_incremental_LiteRate.txt -TBP #innovation linear
#Burn-in set to 20%, meaning the first 20% of sampeles are discarded
python3 plotRJforward.vs.py ~/literate_mcmc_logs/ -TBP -burnin 0.2
Note: When multiple LiteRate files are plotted together, the R object names become overwritten. It may be necessary to re-name these objects manually in order to plot in R (as was done here).
setwd("~/Dropbox/cult_tempo/") #adjust to GitHub when ready
library(scales)
load("neutural_pop_data.RData") #Empirical Replicate Data
source("./neutral_population/neutral_population_literate_raw.R") #Modified (Re-named) Raw Data from LiteRate
#Population Growth - Linear
exp_richness_linear<-apply(matrix(c(pop_sizes$linear,mu),nrow=ts,ncol=2),1,function(x){return(exp.std.div(x[1],x[2]))})
exp_orig_rate_linear <- pop_sizes$linear*mu / exp_richness_linear
exp_ex_rate_linear <- pop_sizes$linear*mu / exp_richness_linear
exp_net_rate_linear <- exp_orig_rate_linear - exp_ex_rate_linear
exp_longevity_linear <- 1 / exp_ex_rate_linear
exp_turnover_linear <- exp_orig_rate_linear / exp_ex_rate_linear
#Population Growth - Logistic
exp_richness_logistic<-apply(matrix(c(pop_sizes$logistic,mu),nrow=ts,ncol=2),1,function(x){return(exp.std.div(x[1],x[2]))})
exp_orig_rate_logistic <- pop_sizes$logistic*mu / exp_richness_logistic
exp_ex_rate_logistic <- pop_sizes$logistic*mu / exp_richness_logistic
exp_net_rate_logistic <- exp_orig_rate_logistic - exp_ex_rate_logistic
exp_longevity_logistic <- 1 / exp_ex_rate_logistic
exp_turnover_logistic <- exp_orig_rate_logistic / exp_ex_rate_logistic
#Population Growth - Exponential
exp_richness_exponential<-apply(matrix(c(pop_sizes$exponetial,mu),nrow=ts,ncol=2),1,function(x){return(exp.std.div(x[1],x[2]))})
exp_orig_rate_exponential <- (pop_sizes$exponetial*mu) / exp_richness_exponential
exp_ex_rate_exponential <- (pop_sizes$exponetial*mu) / exp_richness_exponential
exp_net_rate_exponential <- exp_orig_rate_exponential - exp_ex_rate_exponential
exp_longevity_exponential <- 1 / exp_ex_rate_exponential
exp_turnover_exponential <- exp_orig_rate_exponential / exp_ex_rate_exponential
#Population Growth - Cycle
exp_richness_cycle<-apply(matrix(c(pop_sizes$cycle,mu),nrow=ts,ncol=2),1,function(x){return(exp.std.div(x[1],x[2]))})
exp_orig_rate_cycle <- N*mu / exp_richness_cycle
exp_ex_rate_cycle <- N*mu / exp_richness_cycle
exp_net_rate_cycle <- exp_orig_rate_cycle - exp_ex_rate_cycle
exp_longevity_cycle <- 1 / exp_ex_rate_cycle
exp_turnover_cycle <- exp_orig_rate_cycle / exp_ex_rate_cycle
#TBD
#TBD